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Abstract 

Gaussian wave packets (GWPs) are well suited as basis functions to describe the 
time evolution of arbitrary wave functions in systems with nonsingular smooth 
potentials. They are less so in atomic systems on account of the singular behavior 
of the Coulomb potential. We present a time-dependent variational method that 
makes the use of GWPs possible in the description of propagation of quantum states 
also in these systems. We do so by a regularization of the Coulomb potential, and by 
introduction of a fictitious time coordinate in which the evolution of an initial state 
can be calculated exactly and analytically for a pure Coulomb potential. Therefore 
in perturbed atomic systems variational approximations only arise from those parts 
of the potentials which deviate from the Coulomb potential. The method is applied 
to the hydrogen atom in external magnetic and electric fields. It can be adapted to 
systems with definite symmetries, and thus allows for a wide range of applications. 



1 Introduction 



Wave packets in atomic systems can be excited experimentally, e.g., with mi- 
crowaves [T|2|3] or short laser pulses p|5j . Theoretically, the time evolution 
of wave packets can be calculated accurately by numerically solving the time- 
dependent Schrodinger equation [6], or by using approximation methods, such 
as semiclassical [7] or time-dependent variational ^ methods. The topic of 
wave packet dynamics in systems with Coulomb interactions covers a large 
body of problems ranging from atomic physics to physics of solid state, where 
Coulomb interaction plays an important, often crucial, role. In many-body 
physics, in particular, in the physics of solid state, theoretical methods well 
suited for studying the effects stemming from Coulomb interactions are still 
lacking. The majority of the available methods, e.g., the method of pseudopo- 
tentials in atomic physics and the Fermi and the Luttinger liquid theories for 
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solid conductors, are basically indirect and substantiated neither from the the- 
oretical nor from the experimental side. For this reason they still remain, to a 
certain extent, disputable. The time- dependent variational principle (TDVP) 
applied to Gaussian wave packets (GWPs) lead to exact results for the har- 
monic oscillator potential. GWPs have turned out to be also well suited for 
describing the time evolution of arbitrary wave functions in smooth and nearly 
harmonic potentials piTO] but they are bound to fail for atomic systems be- 
cause of the singularity of the Coulomb potential. It is the objective of this 
Letter to make the GWP method applicable to the description of the time 
evolution of arbitrary quantum states also in these systems. 

For the one- dimensional (ID) Coulomb potential, attempts already have been 
made pT|12|13] to use GWPs as trial wave functions, based on a local harmonic 
approximation. For the full 3D Coulomb potential, which we will consider, a 
way to remove the singularity is, as is well known, the transformation to 4D 
Kustaanheimo-Stiefel (KS) coordinates [TillTSj . which converts the Coulomb 
potential into a sum of two 2D harmonic oscillator potentials, adapted to 
the use of GWPs, but also introduces an additional constraint on the wave 
functions. The regularization implies a fictitious time coordinate. Various ap- 
proaches have been made to construct coherent states for the hydrogen atom 
jl6fl7fl8|19|20j in the fictitious time in analogy with the coherent states of 
the harmonic oscillator. These approaches construct the coherent states as the 
eigenstates of the lowering operators associated with the harmonic potential. 

We will present a variant of the GWP method in coordinate space, which de- 
scribes wave packet propagation in the Coulomb problem exactly. Therefore 
it is only deviations from the Coulomb potential which require a variational 
treatment. As prime examples we will apply the method to wave packet prop- 
agation in the hydrogen atom in a magnetic field, and in crossed electric 
and magnetic fields. Both systems have attracted considerable attention over 
the past decades because classically they exhibit a transition from regular to 
chaotic motion and thus can be used in the search for quantum signatures of 
chaos pni^ . If the method can be further extended to larger systems with 
more degrees of freedom it will allow for a wide range of future applications 
in different branches of physics. 



2 Restricted Gaussian wave packets 



The Hamiltonian for an electron under the combined action of the Coulomb 
potential and external perpendicularly crossed electric and magnetic fields has 
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the form (in atomic units, with Fq = 5.14 x 10^ V/cm, Bq = 2.35 x 10^ T) 

^^3 = ^p' - ^ + ^Bh + Ib^x' + y') + Fx . (1) 

Here we have assumed that the electric field is oriented along the x axis, 
and the magnetic field along the z axis. We regularize the singularity of the 
Coulomb potential by switching to KS coordinates u with x = uiUs — M2M4, 
y = uiUi + U2U3, and z = ^{ul + u1 — u\ — M4). Introducing scaled coordinates 

and momenta u "^cff^u, ^cs'^Vu one obtains 



V^ = 2neff^, (2) 



where the scaled potential V depends on the parameters 

a = -nl^E, l3 = nlf,B, C = nl^F , (3) 

which can be chosen constant. Eq. ([2]) is an eigenvalue problem for the effec- 
tive quantum number Ues, and for any quantized riefr the energy E and field 
strengths B and F of the physical state are obtained from Eq. Q. In KS 
coordinates physical wave functions must fulfill the constraint 

{U2P1 - U1P2 - U4P3 + U3P4) iJ = . (4) 

For a = 1/2 and vanishing external fields {P = 0, ( = 0) Eq. ([2]) describes 
the 4D harmonic oscillator with V = ^u^, and n^s = n becomes the principal 
quantum number of the field-free hydrogen atom. 



Eq. ([2]) can be extended to the time-dependent Schrodinger equation in a 
fictitious time r by the replacement 2nes — > i^, viz. 

= (2?' + V)iJ = Hi;. (5) 

When the TDVP is used to solve Eq. ([5]) the wave function ip depends on a set 
of appropriately chosen parameters whose time-dependences are obtained by 
solving ordinary differential equations. As the regularized Hamiltonian with- 
out external fields ([2]) becomes that of a harmonic oscillator basis trial wave 
functions in the form of GWPs 

g{A, q, TT, 7) = e^[("-q)^("-q)+-("-q)+7] (q) 

with time-dependent parameters are a natural choice. In ([6]) A designates a 
complex symmetric 4x4 width matrix with positive definite imaginary part, 
TT and q are the expectation values of the momentum and position opera- 
tor, and the phase and normalization are given by the complex scalar 7. In 
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KS coordinates physical wave functions must fulfill the constraint (jl]). Insert- 
ing the ansatz into (jl]) leads to restrictions for the admissible variational 
parameters, viz. q = 0, tt = 0, and the special form of the width matrix 



A 



\ 



ttx CLf 

ay -a, 
\^ tty —a^ a,^ J 



which depends only on four parameters (a 



(7) 



,) ^M^. The "restricted 



Gaussian wave packets" obeying Eq. (jl]) are located around the origin with 
zero mean velocity and thus at first glance might not appear appropriate for 
dynamical calculations. However, they are the key for both the exact ana- 
lytical derivation of the fictitious time wave packet dynamics in the field-free 
hydrogen atom and the time-dependent variational approach to the perturbed 
atom. 

The restricted GWPs are not a complete basis set for the four- dimensional 
harmonic oscillator but they are in the 3D space, i.e., any physically allowed 
state can be expanded in that basis. This can be verified by transforming the 
restricted GWP in KS coordinates back into 3D Cartesian coordinates, 



9{y) 



3i(uAu+7) 

^i[{afj,+av)r-\-{a^-av)z+2axX-\-2ayy-\-^] 
3i(prr+p-x+7) 



^b) 

3c) 



In (IHcl) the set of parameters {afj,,au,ax,ay,'j) is replaced with an equivalent 
set y = (pr,P,7) with 



Pr 



a^^ + a^, p = {px, py, p^) = (2a^, 2ay, - a^) . 



(9) 



For Pr = and real valued parameters Px,Py,Pz the restricted GWP in Carte- 
sian coordinates flScj) reduces to a plane wave. Since plane waves form a com- 
plete basis we have the result that the restricted GWPs (jS]) are also complete, 
or even over-complete. 



3 Time-dependent variational principle 



The propagation of the wave packets is investigated by applying the TDVP. 
Briefly, the TDVP of McLachlan ^25j, or equivalently the minimum error 
method [2S], requires to minimize the deviation between the right-hand and 
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the left-hand side of the time-dependent Schrodinger equation with the trial 
function inserted. The quantity 



I = \\i(l){T) - Hij{T)\\^ = min (10) 

is to be varied with respect to (p only, and then ip = (/) is chosen, i.e., for any 
time r the fixed wave function iP{t) is supposed to be given and its admissible 
time derivative iplr) is determined by the requirement to minimize I. As trial 
functions we consider superpositions of N restricted GWPs ( l8a|l . i.e., 

N N 

V;(r)=^(z(r)) = 5:^(y^) = E/, (11) 

k=l k=l 

which are parameterized by a set of 5A^ time-dependent complex parameters 
z = (y'^, k = 1, . . . , A^) [instead of 15 complex parameters when using the 
most general superposition of Gaussian wave packets ([6]) in 4D coordinate 
space]. The equations of motion for the variational parameters z(r) are ob- 
tained as 

itiA'^-v^, (12b) 



.k 
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where we have introduced the time-dependent scalars Vq and matrices V2, 
which induce couplings between the restricted GWPs. Since the special struc- 
ture of the matrices in Eq. is maintained in the squared matrices 
(A^)^, that structure carries over to the 4x4 complex symmetric matrices 
V2. Therefore, they have only four independent coefficients (V^, V^, V^, Vy), 
in the notation of Eq. ([7]). The parameters and V2 are calculated at each 
time step by solving a 5A^ dimensional set of linear equations. All integrals 
required for the setup of that linear system have the form (5''|/(u, p^)!^''^), 
with /(u, p„) a polynomial in the KS coordinates and momenta, and can be 
calculated analytically [21] . 

For the field-free hydrogen atom one finds f g = and V2 = 1, i.e., the 
equations of motion ([T^ simplify to the uncoupled equations A = —2 A? — |l, 
7 = z tr y4 for the parameters of each each basis state. These equations can be 
solved analytically [23] and yield for the time evolution of the restricted GWP 
the explicit form 



with 



Z(r) = p°-x + p°rcos2r + ^[(p°)2-(p°)2-l]sin2r, (14a) 
M{t) = l + [(j9°)2-(p0)2](l-cos2r)+p°sin2r , (14b) 
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and where and p° are the parameters ([9]) of the initial GWP at time r = 0. 
This is an important result for the field-free hydrogen atom: The time evolution 
of a restricted GWP (IScI) can be calculated analytically, and takes the compact 
form fll3p . which is a periodic function of the fictitious time r with period vr. 
In the physical time wave packets disperse in the hydrogen atom. By contrast, 
the wave packets in the fictitious time show an oscillating behavior, with no 
long-time dispersion in r. 



4 Gaussian wave packet dynamics 

We now investigate the propagation of 3D Gaussian wave packets which are lo- 
calized around a given point xq with width a in coordinate space, and around 
Po in momentum space. As mentioned above any physical state can be ex- 
pressed in terms of the complete basis set of the restricted GWPs The 
Fourier decomposition of the initial 3D GWP has the form 

?/' = (27r(T2)-^/^exp 

/ 2 \ 3/4 

= (l^j / ^'pe-'^'(p-p»)'-p-»^7(y) , (15) 

where the g{y) are the restricted GWPs ( IHcl) for the set of parameters y given 
as {pr = 0,p,7 = 0). 

In numerical computations it is convenient to approximate the initial 3D Gaus- 
sian wave packet by a finite number of restricted GWPs rather than using the 
integral representation f[T^ . This is most efficiently achieved by evaluating the 
integral in f[T^ by a Monte Carlo method using importance sampling of the 
momenta. The initial wave packet then reads 

1 ^ 

^ = (2vra2)-3/4 ^ g (yfc) , (16) 

k=i 

with y'^ = (ie, p*^ — 2exo/|xo|, 0), and the p^ distributed randomly according to 
the normalized Gaussian weight function w(p) = (ct^/tt)^/^ exp{— cr^(p— po)^}. 
A small e > has been introduced for damping of the restricted GWPs at large 
radii r, which is convenient in numerical computations. The wave function ip 
in Eq. flTBl) is an approximation to the 3D Gaussian wave packet ffTB]) . and 
the accuracy depends on how many restricted GWPs are included. However, 
it is important to note that a localized wave packet can be described even 
with a rather low number N of restricted GWPs. The time propagation of 
an initial state (fTBjl in the fictitious time r is now obtained exactly and fully 



r (x-xo)^ 1 

I i^+^Po-(^-^o)| 
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Fig. 1. Spectra with (a) even and (b) odd z parity of the Hamiltonian (l2|) with 
a = 0.5, P = 0.05, C = 0-01 obtained from the propagation of two different 3D 
GWPs. Green and red hne: xq = (6, 0, 0), po = (0, ±1/^2, l/V^), respectively. The 
eigenvalues are extracted from the autocorrelation function by Fourier transform. 
The peak positions agree very well with the numerically exact eigenvalues of the 
effective quantum number marked by blue lines. 

analytically by replacing the initial restricted GWPs g{'y'') in Eq. ( |T6i) with 
the corresponding time-dependent solutions (fT3|) . Results for the wave packet 
propagation in the field-free hydrogen atom are given elsewhere |23j . 

Here we present example calculations for the hydrogen atom in external fields. 
In crossed electric and magnetic fields the propagation of 3D GWPs is com- 
puted for the time-dependent Schrodinger equation ([21) with parameters a = 
0.5, P = 0.05, and ( = 0.01 in Eq. ([3]). The choice of an appropriate initial 
state i/j{0) is very important for the successful application of the TDVP. We 
achieved optimal results by choosing a 3D initial Gaussian wave packet in 
physical Cartesian coordinates as given in f|T5l) . The external fields lead to 
couplings between the basis states and the time- dependence of the variational 
parameters must be determined by the numerical integration of Eq. f[T^ . For 
better numerical performance we resort to the TDVP with constraints [27]. 

Once a time-dependent wave packet (ITTl) is determined the eigenvalues ries of 
the stationary Schrodinger equation and thus a quantum spectrum of the 
hydrogen atom in external fields can be obtained by frequency analysis of the 
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Fig. 2. (a) Real part of the even z-parity autocorrelation function and (b) spectrum 
of the diamagnetic hydrogen atom extracted from the signal C^(t) by harmonic 
inversion at the field-free ionization threshold E = 0, (3 = n'^gB = 0.5. For com- 
parison the exact eigenvalues in (b) are marked by blue lines. They are in excellent 
agreement with the variational results (red lines). 

time signal 

C(r) = (^(0)|^(r)) = ^c,e-2^'^o^, (17) 

j 

with the amplitudes Cj depending on the choice of the initial wave packet. 
In perpendicularly crossed fields the z parity is conserved. Spectra with even 
and odd z parity obtained from the Fourier transforms of the autocorrelation 
functions C^(r) = {tp"^ {0)\ip'^ (r)) of the parity projected wave packets are 
shown in Fig. [H The green and red lines result from the propagation of two 
different initial 3D GWPs with a = 3.5, e = 0.15, the same initial position 
but different initial mean momenta. = 41 and iV = 31 basis states were 
coupled in the calculations. The line widths, i.e., the resolution of the spectra 
is determined by the length of the time signal Tmax- The eigenvalues obtained 
by numerically exact diagonalizations of the stationary Hamiltonian ([2]) are 
shown by the blue lines. The line-by-line comparison shows very good agree- 
ment between the exact spectrum and the results obtained from the wave 
packet propagation. The amplitudes of levels indicate the excitation strengths 
of states with higher or lower angular momentum by the two initial wave 
packets rotating clockwise or anticlockwise around the 2;-axis. 

The method presented can be especially adapted to systems with, e.g., cylin- 
drical or spherical symmetries. For the hydrogen atom in a magnetic field we 
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consider the very challenging regime around the field-free ionization threshold 
E = where the Coulomb and the Lorentz force are of comparable strength, 
resulting in a fully chaotic classical dynamics. The real part of the even z- 
parity autocorrelation function is shown in Fig. [2](a). A number of = 90 
basis states was used in the computation. The eigenvalues are extracted from 
the signal C+(r) by the high-resolution harmonic inversion method [28] and 
drawn in Fig. [2]^b). The agreement between the eigenvalues computed varia- 
tionally (red lines) and the numerically exact results (blue lines) is very good. 
Some lines are lacking in the variational computation because of a nearly zero 
overlap of the respective eigenstates Irigfr) and the initial GWP, but can be 
revealed by choosing different initial GWPs. 



5 Conclusion 

In this Letter we have extended the Gaussian wave packet method in such way 
that it can also be applied to quantum systems with singular Coulomb po- 
tentials. We have shown that the evolution in fitctious time can be calculated 
analytically in the pure quantum Coulomb problem. Therefore in applying 
the time-dependent variational principle to the description of time evolution 
of wave packets in perturbed atomic systems approximations arise only from 
the non-Coulombic parts of the potentials. The method can be adapted to 
special symmetries, such as axisymmetric or spherical, and opens the way to 
a wide range of applications in systems with Coulomb potentials. 
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